function dxx = fun_crtbp_stm_(t , xx0 , auxdata)
mu = auxdata.mu;
mu1 = 1 - mu ;
mu2 = mu ;

x = xx0(1 : 6); % state vector
Phi = reshape(xx0(7 : 42) , 6 , 6); % stm

r2 = (x(1) + mu2)^2 + x(2)^2 + x(3)^2; % r: distance to m1, LARGER MASS
R2 = (x(1) - mu1)^2 + x(2)^2 + x(3)^2; % R: distance to m2, smaller mass

r3 = r2^1.5;
r5 = r2^2.5;
R3 = R2^1.5;
R5 = R2^2.5;

% effective potential U(x,y,z)
Ux = - x(1) + mu1 * (x(1) + mu2) / r3 + mu2 * (x(1) - mu1) / R3 ;
Uy = - x(2) + mu1 * x(2) / r3 + mu2 * x(2) / R3 ;
Uz = mu1 * x(3) / r3 + mu2 * x(3) / R3 ; % ADD

Uxx = -1 + (mu1 / r3) * (1 - (3 * (x(1) + mu2)^2 / r2)) + (mu2 / R3) * (1 - (3 * (x(1) - mu1)^2 / R2));
Uyy = -1 + (mu1 / r3) * (1 - (3 * x(2)^2 / r2)) + (mu2 / R3) * (1 - (3 * x(2)^2 / R2));
Uzz = (mu1 / r3) * (1 - (3 * x(3)^2 / r2))+(mu2 / R3) * (1 - (3 * x(3)^2 / R2));

Uxy = -(mu1 / r5) * 3 * x(2) * (x(1) + mu2) - (mu2 / R5) * 3 * x(2) * (x(1) - mu1);
Uxz = -(mu1 / r5) * 3 * x(3) * (x(1) + mu2) - (mu2 / R5) * 3 * x(3) * (x(1) - mu1);
Uyz = -3 * mu1 * x(2) * x(3) / r5 - 3 * mu2 * x(2) * x(3) / R5;

% jacobian of dynamics
Df = [0 , 0 , 0 , 1 , 0 , 0;
    0 , 0 , 0 , 0 , 1 , 0;
    0 , 0 , 0 , 0 , 0 , 1;
    -Uxx , -Uxy , -Uxz , 0 , 2 , 0;
    -Uxy , -Uyy , -Uyz , -2 , 0 , 0;
    -Uxz , -Uyz , -Uzz , 0 , 0 , 0];
dPhi = Df * Phi;

dxx = zeros(42 , 1);
dxx(1 : 6) = [x(4);
    x(5);
    x(6);
    2*x(5) - Ux ;
    -2*x(4) - Uy ;
    -Uz];
dxx(7 : 42) = dPhi(:);

end
% -------------------------------------------------------------------------